% This subroutine generates ability distribution and labor mobility
% barriers.
%
% Used in the paper "Land Insecurity and Mobility Frictions," by Tasso
% Adamopoulos, Loren Brandt, Chaoran Chen, Diego Restuccia, and Xiaoyun
% Wei, prepared for publication at the Quarterly Journal of Economics.
%
% Last modified: March 3, 2024


sH    = exp(norminv(HA(:,4),-P.sigma_sH^2/2,P.sigma_sH));
hH    = exp(norminv(HA(:,5),-P.sigma_hH^2/2,P.sigma_hH));

ind_s   = ones(P.N_sim,P.N_ind);
ind_h   = ones(P.N_sim,P.N_ind);

ind_xir = zeros(P.N_sim,P.N_ind);
ind_xiu = zeros(P.N_sim,P.N_ind);

for i = 1:P.N_ind
    omega = exp(norminv(HA(:,8+(i-1)*5),-P.sigma_omega^2/2,P.sigma_omega)); 
    
    sI    = exp(norminv(HA(:,6+(i-1)*5),-P.sigma_sI^2/2,P.sigma_sI));
    hI    = exp(norminv(HA(:,7+(i-1)*5),-P.sigma_hI^2/2,P.sigma_hI));
    
    h     = omega.*hH.*hI;
    s     = omega.*sH.*sI.*h.^P.omega;
    
    s     = s.*exp(P.mu_sI_y.*(1-dummy_old(:,i)));
    h     = h.*exp(P.mu_hI_y.*(1-dummy_old(:,i)));
    
    ind_s(:,i)   = s;
    ind_h(:,i)   = h;
    
    ind_xir(:,i) = exp(P.zeta_xir*log(ind_h(:,i))+...
        norminv(HA(:,9+(i-1)*5),P.mu_xir_o-P.sigma_xir^2/2,P.sigma_xir))./...
        (1+exp(P.zeta_xir*log(ind_h(:,i))+...
        norminv(HA(:,9+(i-1)*5),P.mu_xir_o-P.sigma_xir^2/2,P.sigma_xir)))...
        .*dummy_old(:,i) + ...
        exp(P.zeta_xir*log(ind_h(:,i))+...
        norminv(HA(:,9+(i-1)*5),P.mu_xir_y-P.sigma_xir^2/2,P.sigma_xir))./...
        (1+exp(P.zeta_xir*log(ind_h(:,i))+...
        norminv(HA(:,9+(i-1)*5),P.mu_xir_y-P.sigma_xir^2/2,P.sigma_xir)))...
        .*(1-dummy_old(:,i));
    
    ind_xiu(:,i) = exp(P.zeta_xiu*log(ind_h(:,i))+...
        norminv(HA(:,10+(i-1)*5),P.mu_xiu_o-P.sigma_xiu^2/2,P.sigma_xiu))./...
        (1+exp(P.zeta_xiu*log(ind_h(:,i))+...
        norminv(HA(:,10+(i-1)*5),P.mu_xiu_o-P.sigma_xiu^2/2,P.sigma_xiu)))...
        .*dummy_old(:,i) + ...
        exp(P.zeta_xiu*log(ind_h(:,i))+...
        norminv(HA(:,10+(i-1)*5),P.mu_xiu_y-P.sigma_xiu^2/2,P.sigma_xiu))./...
        (1+exp(P.zeta_xiu*log(ind_h(:,i))+...
        norminv(HA(:,10+(i-1)*5),P.mu_xiu_y-P.sigma_xiu^2/2,P.sigma_xiu)))...
        .*(1-dummy_old(:,i));
end

 